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Abstract 

The ALICE and CMS data on the multiphcity distributions are compared with 
the lower energy data and with the results from the 8.142 version of the PYTHIA 
MC event generator with two tunings. The ALICE data for moments are used 
to calculate the factorial cumulants. It is suggested that the data on moments or 
cumulants are well suited to specify the optimal tuning of the model parameters. 
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1 Introduction 



With the advent of LHC data it became possible to investigate the multiphcity distri- 
butions at the CM energies beyond 2 TeV. It is very interesting to check how well the 
default versions of MC generators describe the minimum bias events in this energy range. 
In particular, it is tempting to look for a best tuning of the model parameters using only 
the data from the multiplicity distributions. 

In a recent note [Ij we have discussed the energy dependence of the central density 
(defined by the average charged multiplicity in a central bin in pseudorapidity). We have 
shown that, contrary to some claims, the fast increase with energy observed in the ALICE 
data is not unexpected. In fact, we found that the default version of PYTHIA 8.135 
generator [3], jl] predicts a too fast increase, but with some tuning the data may be 
well described. The moments for three selected pseudorapidity bins were also compared 
with data [5]. The qualitative agreement was observed and the spread of results for two 
different tunings was surprisingly small. 

In this paper we discuss the data on the multiplicity distributions listed above as well 
as the data from the CMS experiment [6] compared with some lower energy data [7], |8], 
[9]. We use here the new version of PYTHIA 8.142 [10] with the default tuning and with 
the tuning used by default in the older version 8.135. Let us note that the PYTHIA 
8.142 version is significantly changed compared to its earlier versions (PYTHIA 8.107 and 
8.135) which we have used in a previous publication pj. The changes, concerning mainly 
the final state radiation, are motivated by the discrepancy between the Tevatron data for 
"underlying event" and the model results. The qualitative conclusions of our note [1] are, 
however, unchanged. 

In the next section, we compile for convenience the formulae defining the moments, 
factorial moments and factorial cumulants of the multiplicity distributions. Then we recall 
the results from Ref. [1] and compare them with the factorial cumulants calculated from 
the data and MC generators. We will see that that the scaled factorial cumulants Kq, 
calculated from the published values of the standard scaled moments Cq and the average 
multiplicities rT, exhibit much smaller spread than suggested by the published values of 
the uncertainties of Cq and n. Thus it would be more reasonable to use the scaled factorial 
cumulants of the multiplicity distributions to test the specific models of the high energy 
collisions. Finally, we discuss the energy dependence of the average multiplicity and 
second moment in a wider energy range, as presented in the CMS paper [6]. The last 
section contains some conclusions and the outlook. 



2 Moments and cumulants 

The multiplicity distributions are often parametrized in terms of moments. This facilitates 
the comparison with models and allows for a simple description of the energy dependence. 
The crucial problem is the proper choice of the set of moments to be used. A standard 
first choice is to use simple power moments defined by 

= T.n''P{n) 

or their scaled version 

_ nt 

These moments are easy to calculate and (for moderate q) depend quite uniformly 
on the probabilities, although obviously the lowest multiplicities are suppressed, and the 



high multiphcity tail is enhanced. The use of n and a few lowest Cq moments allows to 
parametrize the multiplicity distribution quite satisfactorily. In the high energy (high n) 
limit these moments allow to describe the "KNO scaling function" \E'(x) [11] defined by 

\E'(x) = Jim P{n)/n 

n—^oo 

where x = n/n. Obviously, in this limit 

However, for the distributions in small phase space bins the power in the denominator 
of the formula for the scaled moments results in the development of trivial singularities. 
Thus since some time one uses more often the factorial moments. If we define the factorial 
quotient 

n! 



" {n-q)\ 

the corresponding standard and scaled factorial moments are, respectively 



and 



Uq = T,nqP(n) 



The factorial moments of the order q are the integrals of the q-particle densities for 
identical particles. Obviously, nT = n is the integral of the single particle density. For the 
smooth phase space distributions the scaled factorial moments behave smoothly for the 
bin size decreasing to zero, and the possible power increase is a signal for intermittency 
|12) . However, the drawback of the definition of the higher factorial moments is their 
independence on the lower end of the multiplicity distribution. Moreover, the scaled 
factorial moments of different order are strongly correlated. 

Therefore it is preferable to parametrize the multiplicity distributions by the factorial 
cumulants. They are defined in a compact way by the generating function 

G{z) = S2"P(n). 

The factorial cumulants fq (called also "Mueller coefficients" [13]) are defined by 

d'^{lnG{z)) 



z=l 



dz'i 

to be compared with an analogous definition of the factorial moments 

d^G{z) 



^ dzi 



z=l 

The scaled factorial moments are defined in a usual way 



K,--q. 



By definition, for g = 1 all the scaled moments Cg, Fq and Kq are equal one. Another name 
for the factorial cumulants fq is "the correlation integrals", as they are the integrals of 



the correlation functions of the order q. Therefore they measure the genuine muhiparticle 
correlations. For the uncorrelated emission all the factorial cumulants for g > 1 vanish. 
If there are only two-particle correlations, fq = for q > 2. The values of the factorial 
cumulants of different orders are uncorrelated. 

One may add that another set of moments was advocated [T^ 

Hq = KjFq = fjn-q. 

These moments were shown to have strongly reduced statistical uncertainties even for the 
order up to g = 10. However, for g < 5 they are not very practical to use, as their values 
quickly decrease with increasing q. 

A practical difficulty in using the factorial cumulants is the complexity of the formulae 
for their errors, or, more precisely, for their statistical uncertainties. For the average 
multiplicity one uses a simple estimate of the uncertainty 

Arr = L)/ViV 

where D is the dispersion, and is the total number of measured events. Analogous 
simple formulae exist for the higher standard moments n"?. However, for the scaled mo- 
ments, and especially for cumulants, the corresponding formulae are more complicated. 
Moreover, usually they overestimate significantly the observed spread of experimental 
results. 

There is a simple explanation of this fact. The formula for An was obtained from a 
simple prescription for the statistical uncertainty of a parameter 

= ^T.{dA/dNnY ■ (AAr„)2 

with AA^„ = ^/Nn- This prescription results from the assumption that the measured num- 
bers of events with different multiplicities A^„ are uncorrelated, and their errors are purely 
statistical. These assumptions were reasonable e.g. for the hydrogen bubble chamber 
experiments, where the full solid angle was available for the measurements of tracks, and 
the multiplicity of charged particles (always even for charged beam) was unambiguously 
measured. This is certainly not the case for a colliding beam experiment with electronic 
detectors, where the multiplicity distribution is measured in the restricted bin of phase 
space. The uncertainty of a measurement of the variables defining this bin as well as the 
effects of the track splitting and joining due to the imperfection of the detector result in 
non-statistical errors and in the strong correlations between the numbers of events with 
different multiplicities. Neither a simple formula for An presented above, nor the compli- 
cated formulae for the uncertainties of the scaled factorial cumulants (derived from the 
same prescription for the statistical uncertainties of the numbers of events) are reliable. 

Therefore the realistic estimate of the uncertainties of the parameters of the mul- 
tiplicity distribution requires the full knowledge of the detector and, in particular, the 
measurement of the correlation matrix for the multiplicities. This can be done only by 
the authors of the experiment. Readers cannot translate them reliably into a different set 
of parameters, since their uncertainties will be unknown. 

3 Moments and factorial cumulants 

In a recent paper [5j the ALICE collaboration has presented the values of average multi- 
plicities and scaled moments Cq for g = 2,3,4 at two CM energies 0.9 and 2.36 TeV for 
three choices of the central pseudorapidity bin widths: A?] < 1, A?] < 2 and A?] < 2.6. 



In Table 1 and in Fig.l we show the experimental values of the n from the "non-single- 
diffractive" (NSD) ALICE data at 900 GeV and 2.36 TeV and the corresponding values 
calculated from the PYTHIA 8.142 default version. For each point we have generated 
10^ events. In all tables the numbers in parentheses denote the statistical and systematic 
errors. 

Table 1: Average multiplicities for three choices of rapidity bins from ALICE and two 
versions of PYTHIA 8.142 at 0.9 and 2.36 TeV. 



rj range 


ALICE 0.9 TeV 


P8.142 


P8.142/135 


ALICE 2.36 TeV 


P8.142 


P8.142/135 


1 7? |< 0.5 


3.60(2)(11) 


3.57 


3.87 


4.47(3)(10) 


4.45 


4.75 


1 V l< 1-0 


7.38(3)(17) 


7.27 


7.88 


9.08(6)(29) 


9.04 


9.66 


hl<l-3 


9.73(12)(19) 


9.57 


10.35 


11.86(22)(45) 


n.89 


12.68 
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Figure 1: The average multiplicity from the ALICE data at 0.9 and 2.36 TeV (asterisks with error 
bars), from PYTHIA 8.142 (open squares and circles) and from PYTHIA 8.142/135 (full squares and 
circles) as a function of the pseudorapidity bin width 

The data for the average multiplicities agree well with the model (remember that we 
are using the default version of PYTHIA 8.142 without any tuning). The agreement is 
significantly better than that for the central densities of charged particles for inelastic 
events with at least one particle in the central bin, measured by ALICE at 0.9, 2.36 and 
7 TeV [2], although in this case the increase with energy is also reasonably well described, 
as seen in Table 2 and in Fig. 2. 

We have repeated the same calculations for the PYTHIA 8.142 with a different set of 
the model parameter values: the default values from the PYTHIA 8.135 version are taken. 
Let us remind here that the tuned parameters refer to the formulae used in the description 
of multiple scattering. The regularization of the (divergent) QCD cross section is done 
by the introduction of a factor 




Table 2: Central density: data and the results for two versions of PYTHIA 8.142 



Energy (TeV) 


ALICE 


PYTHIA 8.142 d 


PYTHIA 8.142/135d 


0.90 


3.81(1)(7) 


3.58 


3.86 


2.36 


4.70(1)(11) 


4.41 


4.69 


7.00 


6.01(1)(20) 


5.80 


6.14 



A 
V 




Figure 2: Central density: data (asterisks with error bars), PYTHIA 8.142 default (squares) and 
PYTHIA 8.142 with 8.135 tuning (circles) as a function of the CM energy. 



where 

^ , / eCmNoW \ ecmPow 

PTo=pTORef( 

V ecmnej / 

and ecmNow is the CM energy in GeV. The default values in PYTHIA 8.135 are 2.0 for the 
pTORef, 1960.0 for ecmRef and 0.16 for ecmPow. In PYTHIA 8.142 the corresponding 
values are 2.15, 1800.0 and 0.24. Moreover, the default version of PYTHIA 8.142 uses 
a simplified profile of the parton density in the impact parameter given by a Gaussian 
curve, whereas the standard earlier versions were using two Gaussians (with two extra 
parameters for the ratios of their slopes and weights). We have found that in this case the 
results for the average multiplicities are reversed: there is a perfect agreement for central 
densities in the "INEL>0" sample (see Table 2) and a slight overestimation for the NSD 
sample. 

The difference between Tables 1 and 2 shows that the ALICE procedures give for the 
PYTHIA events practically the same average multiplicity in the "NSD" and "INEL>0" 
samples, whereas experimentally the second sample has higher average multiplicity, which 
suggests lower contribution from the diffractive events. Remember that for PYTHIA we 
use the same definition of "NSD" and "INEL>0" events as in the data. This means 
that we generate all the classes of events (non-diffractive, single diffractive and double 
diffractive) and then remove the events which do not satisfy the conditions defined in the 
ALICE procedures. 



The formulae listed in the previous section allow to express the scaled factorial cumu- 
lants in terms of the scaled moments and average multiplicity. For the lowest values of q 
we have 



K2 = F2- 1, F2 = C2 - 1/n, K^ = F^- 3F2 + 2, F3 = C3 - Cs/n + 2/rf 

and 

K4 = F4- 4F3 - 3F| + I2F2 - 6, F4 = C4 - Gcs/n + llcs/n^ - 6/n\ 

For the higher values of q it is more practical to use a recurrence formula expressing Kq 
by Fq and the K moments of the lower order 

In Table 3 we show the values of the Cq and Kq moments at 900 GeV, and in Table 4 
the same results at 2.36 TeV. 



Table 3: Scaled moments and factorial cumulants for three choices of rapidity bin from 
ALICE and two versions of PYTHIA 8.142 at 0.9 TcV. 



T] range 


Cg 


ALICE 


P8.142 


P8. 142/ 135 




ALICE 


P8.142 


P8. 142/135 


\ Tj \ < 0.5 


C.2 


1.96(1)(6) 


1.73 


1.85 


K2 


0.68 


0.45 


0.59 


T] \< 1.0 


C2 


1.77(1)(4) 


1.56 


1.70 


K2 


0.63 


0.42 


0.57 


V l< 1-3 


C2 


1.70(3)(7) 


1.51 


1.65 


K2 


0.60 


0.40 


0.55 


\ r] \ < 0.5 


C3 


5.35(6)(31) 


4.16 


4.93 


K3 


0.82 


0.50 


0.85 


r] |< 1.0 


C3 


4.25(3)(20) 


3.29 


4.11 


K3 


0.66 


0.42 


0.78 


V l< 1-3 


C3 


3.91(10)(15) 


3.04 


3.84 


K3 


0.62 


0.38 


0.73 


T] \< 0.5 


C4 


18.3(4)(1.6) 


12.7 


16.8 


Ki 


1.13 


0.70 


1.39 


1 77 |< 1.0 


C4 


12.6(1)(9) 


8.65 


12.6 




0.82 


0.50 


1.24 


1 r] \< 1.3 


C4 


10.9(4)(6) 


7.60 


11.3 


Ki 


0.57 


0.43 


1.10 



Table 4: Scaled moments and factorial cumulants for three choices of rapidity bin from 
ALICE and two versions of PYTHIA 8.142 at ^.5^ TeV. 



r] range 


Cg 


ALICE 


P8.142 


P8. 142/135 


Kg 


ALICE 


P8.142 


P8. 142/135 


?7 < 0.5 


C2 


2.02(1)(4) 


1.75 


1.90 


K2 


0.80 


0.53 


0.69 


r] \< 1.0 


C2 


1.84(1)(6) 


1.61 


1.76 


K2 


0.73 


0.50 


0.65 


\ V \ < 1-3 


C2 


1.79(3)(7) 


1.56 


1.71 


K2 


0.71 


0.48 


0.64 


?7 < 0.5 


ca 


5.76(9)(26) 


4.31 


5.25 


K3 


1.12 


0.64 


1.08 


1 V \< 1-0 


C3 


4.65(6)(30) 


3.56 


4.46 


K3 


0.88 


0.55 


0.98 


1 V l< 1-3 


C3 


4.35(16)(33) 


3.32 


4.23 


K3 


0.79 


0.51 


0.93 


77 < 0.5 


C4 


20.6(6)(1.4) 


13.4 


18.6 


Ki 


1.77 


0.94 


2.00 


1 77 1< 1.0 


C4 


14.3(3)(1.4) 


9.80 


14.4 


Ki 


0.98 


0.71 


1.68 


1 V l< 1-3 


C4 


12.8(7)(1.5) 


8.75 


13.1 


Ki 


0.83 


0.61 


1.51 



These data show a few simple regularities: 

1. The values of the Cg moments increase with the value of q and with energy, but 
decrease with the increasing pseudorapidity bin width. The average multiplicity, as 
expected, increases with the bin width and energy. 



2. The regularities listed above hold for g = 4 even in the cases, when the experimen- 
tal errors given by the authors exceed the differences between the data for different 
energies or different bin widths. This is not so surprising for the bin width depen- 
dence, as the data are here clearly quite strongly correlated. The presence of a 
similar effect in the energy dependence seems to suggest that the systematic errors 
at two energies are also correlated. 

3. The values of Cq in the two tunings differ by 0.15, 0.9 and 5 for g = 2, 3 and 4, 
respectively. In contrast, the values of Kq differ much less for q > 2. The model 
with the "wrong" tuning is compatible with data for q > 2. This suggests that by 
a more refined tuning one may get a reasonable agreement with data not only for 
average multiplicities, but also for higher moments. Using the factorial cumulants 
we find a smaller spread of the values both in the model and in the data. 

4. The moments of the multiplicity distributions are systematically underestimated in 
the default tuning of PYTHIA 8.142. For the Cq moments the difference between 
the data and the model values is around 0.2, 1 and 5 for q = 2, 3 and 4, respectively, 
and increases weakly with the energy. For the Kg moments the trend is the same, 
but the differences for g > 2 are much smaller: only in one case the difference is 
bigger than 0.5. For the 8.135 tuning the situation is much more involved. All the 
values of the moments are now significantly higher. Whereas for Cq they are still 
lower than in the data, the difference is really significant only for q = 2. For the Kq 
the values are below the data for q = 2, and above the data for g = 4. 

It would be highly desirable to calculate reliably the experimental errors of the scaled 
factorial cumulants to see how significant is the difference between the model and data 
seen in Tables 3 and 4. We have checked that the statistical uncertainties of the model 
results for a given set of parameters are negligible: by increasing the statistics by a factor 
of ten we do not change the values from the Tables by more than a few percent. In 
all cases the observed fluctuations for Cq are negligible compared with the experimental 
errors. However, as noted above, tuning the model parameters allows to change the results 
sufficiently to hope for the agreement with data. 

4 Energy dependence from SPS to LHC 

The CMS collaboration has also measured the multiplicity distributions for non-single- 
diffractive (NSD) events at the CM energies of 0.9, 2.36 and 7 TeV [6]. Contrary to the 
most of published results, where the "NSD" events are defined just by giving the trigger 
conditions, the CMS data are extrapolated and corrected to remove single diffraction. 
Since we are unable to repeat this procedure in detail, the precise comparison with MC 
results is beyond our ability. We have to rely on the effectiveness of the CMS procedure 
and to compare their data with the events generated in PYTHIA as non-single-diffractive. 

The CMS paper refers also to the published data from the lower energies: the CERN 
collider data from UA5 [7J and UAl [8j collaborations and the SPS data from the EHS /NA22 
collaboration |9]. Thus it is possible to check if the energy evolution in the wider range is 
reasonably described by the PYTHIA 8 generator, and what are the differences between 
the different tunings in this range. For the sake of transparency, the UAl data are shown 
only in the figure containing the average multiplicities. 

In Fig. 3 we show the data for the average multiplicity in the range | ^7 |< 2.4 (2.5 for 
lower energies) and the PYTHIA 8.142 results with two tunings described above. At the 
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Figure 3: The average multiplicity in the | 77 |< 2.4 range from the NA22 data (triangle), UA5 data 
(open stars), UAl data (open crosses) and CMS data (black dots). PYTHIA 8.142 predictions for NA22, 
UA5 and CMS with default and 8.135 tunings are shown as x-s and bars. 



CERN collider energies the PYTHIA predictions are shown only with the UA5 trigger 
conditions. We see that the MC results agree quite well with the observed trend of the 
data and the values from the two tunings bracket the experimental results. 
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Figure 4: The C2 moment in the | 77 |< 0.5 and | ?/ |< 2.4 range from the NA22, UA5 and CMS data 
(vertical bars). PYTHIA 8.142 predictions with default and 8.135 tunings are shown as x-s and bars, 
respectively. The points for UA5 and CMS data are connected by dotted lines to guide the eye. The open 
points show the PYTHIA predictions for the pure NSD sample at the UA5 and NA22 energies. 



The situation is much more involved for the moments. In Fig.4 we show the second 
scaled moment C2 for two choices of the rapidity range: | r/ |< 0.5 and | 77 |< 2.4 (2.5 
for lower energies). Apart from the data and MC results for two tunings (which coincide 
practically for lower energies) we show here the MC predictions for "true NSD" events. As 



shown in Fig.4, the MC predictions for the NA22 data (defined by their trigger) are quite 
far from the MC results for "true NSD" events. This casts some doubts on the claims of 
selecting the "true NSD" events from the data by triggers and corrections. For UA5 data 
the difference is smaller. Note that the difference between the PYTHIA predictions for 
UA5 and CMS at the same energy of 0.9 TeV results from the different definition of NSD 
events. 

The agreement of the model with data is unsatisfactory: the moments are underes- 
timated at NA22 and UA5 energies for both rj ranges. However, the differences are not 
too big and the energy dependence is qualitatively correct. The most important feature 
of the results is the sudden increase of differences between the results from two tunings 
at highest energies. This suggests that the reliable measurement of the multiplicity dis- 
tribution at 7 TeV should fix the tuning well enough to allow for a significant test of the 
model from other data. 
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Figure 5: The cs moment in the | |< 0.5 and the | rj \< 2.4 range from the NA22, UA5 and CMS 
data (full circles with error bars). PYTHIA 8.142 predictions with default and 8.135 tunings are shown 
as open squares and triangles, respectively. For transparency, the data points for UA5 and CMS and the 
PYTHIA predictions at 0.9 TeV are slightly shifted to lower and higher energy, respectively. 



In Fig.5 we show the third scaled moment for the rapidity ranges | |< 0.5 and 
I |< 2.4 (2.5 for lower energies). We see that the pattern is similar to that of the second 
moment. For the NA22 data the value of C3 is again much higher than from the smooth 
extrapolation of the higher energy data, suggesting that the trigger does not select well the 
NSD events. The agreement of PYTHIA with low energy data is not quite satisfactory, 
but the energy dependence is qualitatively correct. The LHC data are bracketed by two 
versions of PYTHIA 8.142 tunings, which differ strongly at highest energies. 

In Fig.6 the fourth scaled moment is shown for the rapidity ranges | |< 0.5 and 
I ?7 |< 2.4 (2.5 for lower energies), respectively. Again, the pattern is similar. Note that 
the relative experimental uncertainties are almost the same. This makes the discrepancies 
at lower energies less significant. 

In general, the qualitative agreement of PYTHIA with the energy dependence of the 
mutiplicity distributions should be regarded as acceptable. 
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Figure 6: The C4 moment in the | 77 |< 0.5 and the | 77 |< 2.4 range from the NA22, UA5 and CMS 
data (fuh circles with error bars). PYTHIA 8.142 predictions with default and 8.135 tunings are shown 
as open squares and triangles, respectively. For transparency, the data points for UA5 and CMS and the 
PYTHIA predictions at 0.9 TeV are slightly shifted to lower and higher energy, respectively. 



5 Conclusions and outlook 

We have extended our former analysis of the multiphcity distributions at the LHC ener- 
gies [1] using the new version of MC generator (PYTHIA 8.142 with two tunings) and 
calculated the scaled factorial cumulants from the scaled moments and average multiplic- 
ities for the model and ALICE data. We have also investigated the energy dependence of 
the average multiplicity and the three lowest scaled moments in the wide energy range, 
comparing the PYTHIA results with the NA22, UA5 and CMS data. 

We have found that the fast increase of the central density of charged hadrons at LHC 
energies agrees quite well with the model predictions. The use of factorial cumulants 
should facilitate the fixing of the tuning parameters. The energy dependence of the 
scaled moments is qualitatively well described in the wide range covering the SPS, CERN 
collider and LHC energies. 

There is a large difference between the PYTHIA results with two tunings at the highest 
energies. This suggests that the multiplicity distributions from LHC are well suited to 
fix the tuning of MC generators. Other data could be then compared with the model 
predictions. 
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